library(devtools) 
library(haven)
library(dplyr)
library(emmeans)
library(ggplot2)
library(ggpubrgg)
library(ggrepel)
library(ggeffects)
library(ggthemes)
library(lme4)
library(readr)
library(interactions)
library(patchwork)
library(ggridges)
library(jtools)
library(colorspace)
library(modelsummary)
library(interactions)
library(misty)
library(tidyverse)

LGBTQ_rural_id <- read_dta("LGBTQ_rural_id.dta")

dataframe <- LGBTQ_rural_id%>%
  
mutate(ruralid = as_factor(ruralid), ruralstid = as_factor(ruralstid), place = as_factor(place), region = as_factor(region))

dataframe$V202255x [dataframe$V202255x==-7] <- NA
dataframe$V202255x [dataframe$V202255x==-6] <- NA
dataframe$V202255x [dataframe$V202255x==-5] <- NA
dataframe$V202255x [dataframe$V202255x==-2] <- NA

table(dataframe$V202255x)

dataframe$V202255x <- item.reverse(dataframe$V202255x, min = 1, max = 6)

table(dataframe$V202255x)

names(dataframe)[names(dataframe) == 'V202255x'] <- 'antistate'

table(dataframe$antistate)

dataframe$lgbtindex <- (dataframe$lgbtindex - 0) / (5-0)

#descriptives#

dataframe2 <- subset(dataframe, ruralid=="1")

ggplot(dataframe2, aes(x = placeidsalience)) + geom_bar(aes(y = (..count..)/sum(..count..))) + xlab("Importance of rural identity") + theme(plot.title = element_text(face = "plain")) + scale_x_discrete(labels=c("1" = "Not at all important", "2" = "Not very important", "3" = "Moderately important", "4" = "Very important", "5" = "Extremely important"))  + theme_tufte(base_size=11) +
  scale_y_continuous(labels = scales::percent) + ylab("Percent")

ggsave("Thompson 22-0299 Figure 1.tiff", last_plot(), width = 6, height = 4, dpi = 600, path = "Desktop")

#group-based affect#

model1 <- lm(ft_gays ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model1)

marginsgays <- ggpredict(model1, terms = "ruralid")

ggplot2.1 <- ggplot(marginsgays, aes(x, predicted, color= x)) + geom_point() + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 0.8) + theme_tufte(base_size=12) + 
  ylab("Rating") +
  xlab("Identity") + scale_x_discrete(labels=c("0" = "Non-rural", "1" = "Rural")) + labs(subtitle = "Thermometer: Gays and lesbians") + theme(plot.title = element_text(face = "plain")) + scale_color_manual(values = c("0" = "grey45","1" = "grey45")) + theme(legend.position = "none") + ylim(50, 70) 

emmeans1 <- emmeans(model1, "ruralid")
pairs(emmeans1)

marginssubplace1 <- ggpredict(model1, terms = "place")

marginspid1 <- ggpredict(model1, terms = "pid7")

marginsEvangelical1 <- ggpredict(model1, terms = "Evangelical")

model2 <- lm(ft_trans ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model2)

marginstrans <- ggpredict(model2, terms = "ruralid")

ggplot2.2 <- ggplot(marginstrans, aes(x, predicted, color= x)) + geom_point() + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 0.8) + theme_tufte(base_size=12) + 
  ylab("Rating") +
  xlab("Identity") + scale_x_discrete(labels=c("0" = "Non-rural", "1" = "Rural")) + labs(subtitle = "Thermometer: Transgender individuals") + theme(plot.title = element_text(face = "plain")) + scale_color_manual(values = c("0" = "grey45","1" = "grey45")) + theme(legend.position = "none") + ylim(50, 70) 


emmeans2 <- emmeans(model2, "ruralid")
pairs(emmeans2)

marginssubplace2 <- ggpredict(model2, terms = "place")

marginspid2 <- ggpredict(model2, terms = "pid7")

marginsEvangelical2 <- ggpredict(model2, terms = "Evangelical")

modelsgroupaffect <- list()
modelsgroupaffect [['Gays and lesbians']] <- lm(ft_gays ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
modelsgroupaffect [['Transgender individuals']] <- lm(ft_trans ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)

msummary(modelsgroupaffect, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

Figure2 <- (ggplot2.1 + ggplot2.2)

ggsave("Thompson 22-0299 Figure 2.tiff", last_plot(), width = 8, height = 4, dpi = 600, path = "Desktop")


#LGBTQ rights#

model3 <- lm(lgbtindex ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model3)

marginsrights <- ggpredict(model3, terms = "ruralid")

ggplot(marginsrights, aes(x, predicted, color= x)) + geom_point() + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 0.8) + theme_tufte(base_size=12) + 
  ylab("Index value") +
  xlab("Identity") + scale_x_discrete(labels=c("0" = "Non-rural", "1" = "Rural")) + theme(plot.title = element_text(face = "plain")) + scale_color_manual(values = c("0" = "grey45","1" = "grey45")) + theme(legend.position = "none") + ylim(0.6, 0.75)

ggsave("Thompson 22-0299 Figure 3.tiff", last_plot(), width = 4, height = 4, dpi = 600, path = "Desktop")

msummary(model3, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

emmeans3 <- emmeans(model3, "ruralid")
pairs(emmeans3)

marginssubplace3 <- ggpredict(model3, terms = "place")

#rural ID moderated by place-based identity salience#
#Group based affect#
model4 <- lm(ft_gays ~ ruralid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model4)

marginsgaysmod <- ggpredict(model4, c("placeidsalience", "ruralid"))

ggplot4.1 <- ggplot(marginsgaysmod, aes(x, predicted, color= group)) + geom_point() + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 0.8) + theme_tufte(base_size=12) +
  ylab("Rating") +
  xlab("Importance of place-based identity") + labs(subtitle = "Thermometer: Gays and lesbians") + theme(plot.title = element_text(face = "plain")) + scale_x_continuous(labels=c("1" = "Not at all important", "2" = "Not very important", "3" = "Moderately important", "4" = "Very important",     "5" = "Extremely important")) + theme(legend.title=element_blank()) + geom_smooth(method = "lm") + scale_color_manual(values = c("0" = "grey15","1" = "grey55"), labels = c("0" = "Non-rural", "1" = "Rural")) + theme(axis.text.x = element_text(size = 8)) + ylim(50, 70) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

model5 <- lm(ft_trans ~ ruralid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model5)

marginstransmod <- ggpredict(model5, c("placeidsalience", "ruralid"))

ggplot4.2 <- ggplot(marginstransmod, aes(x, predicted, color= group)) + geom_point() + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 0.8) + theme_tufte(base_size=12) +
  ylab("Rating") +
  xlab("Importance of place-based identity") + labs(subtitle = "Thermometer: Transgender individuals") + theme(plot.title = element_text(face = "plain")) + scale_x_continuous(labels=c("1" = "Not at all important", "2" = "Not very important", "3" = "Moderately important", "4" = "Very important",     "5" = "Extremely important")) + theme(legend.title=element_blank()) + geom_smooth(method = "lm") + scale_color_manual(values = c("0" = "grey15","1" = "grey55"), labels = c("0" = "Non-rural", "1" = "Rural")) + theme(axis.text.x = element_text(size = 8)) + ylim(50, 70) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

modelsgroupaffectintst <- list()
modelsgroupaffectintst [['Gays and lesbians']] <- lm(ft_gays ~ ruralid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
modelsgroupaffectintst [['Transgender individuals']] <- lm(ft_trans ~ ruralid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)

msummary(modelsgroupaffectintst, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

Figure4 <- (ggplot4.1 + ggplot4.2)

ggsave("Thompson 22-0299 Figure 4.tiff", last_plot(), width = 8, height = 4, dpi = 600, path = "Desktop")


#LGBTQ rights#
model6 <- lm(lgbtindex ~ ruralid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model6)

marginsrights <- ggpredict(model6, c("placeidsalience", "ruralid"))

ggplot(marginsrights, aes(x, predicted, color= group)) + geom_point() + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 0.8) + theme_tufte(base_size=12) +
  ylab("Index value") +
  xlab("Importance of place-based identity") + theme(plot.title = element_text(face = "plain")) + scale_x_continuous(labels=c("1" = "Not at all important", "2" = "Not very important", "3" = "Moderately important", "4" = "Very important",     "5" = "Extremely important")) + theme(legend.title=element_blank()) + geom_smooth(method = "lm") + scale_color_manual(values = c("0" = "grey15","1" = "grey55"), labels = c("0" = "Non-rural", "1" = "Rural")) + theme(axis.text.x = element_text(size = 8)) + ylim(0.6, 0.75) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

ggsave("Thompson 22-0299 Figure 5.tiff", last_plot(), width = 4, height = 4, dpi = 600, path = "Desktop")

msummary(model6, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))


#model place id as ordinal#
#Group based affect#
model7 <- lm(ft_gays ~ placeid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model7)

model8 <- lm(ft_trans ~ placeid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model8)

modelsgroupaffectord <- list()
modelsgroupaffectord [['Gays and lesbians']] <- lm(ft_gays ~ placeid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
modelsgroupaffectord [['Transgender individuals']] <- lm(ft_trans ~ placeid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)

msummary(modelsgroupaffectord, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

#LGBT rights#
model9 <- lm(lgbtindex ~ placeid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model9)

msummary(model9, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

#rural ID moderated by place-based identity salience#
#Group based affect#
model10 <- lm(ft_gays ~ placeid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model10)

model11 <- lm(ft_trans ~ placeid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model11)

modelsgroupaffectordint <- list()
modelsgroupaffectordint [['Gays and lesbians']] <- lm(ft_gays ~ placeid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
modelsgroupaffectordint [['Transgender individuals']] <- lm(ft_trans ~ placeid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)

msummary(modelsgroupaffectordint, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

#LGBT rights#
model12 <- lm(lgbtindex ~ placeid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model12)

msummary(model12, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

#individual items#
model13 <- glm(services ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
summ(model13)

model14 <- glm(bathrooms ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
summ(model14)

model15 <- glm(discrimlaws ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
summ(model15)

model16 <- glm(adoption ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
summ(model16)

model17 <- glm(marriage ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
summ(model17)

modelsindividualitems <- list()
modelsindividualitems [['Services']] <- glm(services ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
modelsindividualitems [['Bathrooms']] <- glm(bathrooms ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
modelsindividualitems [['Discrimination laws']] <- glm(discrimlaws ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
modelsindividualitems [['Adoption']] <- glm(adoption ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)
modelsindividualitems [['Same-sex marriage']] <- glm(marriage ~ ruralid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, family = binomial(link = "probit"), data = dataframe, weight = weight_post)

msummary(modelsindividualitems, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

#collapsing rural and small town#
#Group based affect#
model18 <- lm(ft_gays ~ ruralstid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model18)

model19 <- lm(ft_trans ~ ruralstid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model19)

modelsgroupaffectord <- list()
modelsgroupaffectord [['Gays and lesbians']] <- lm(ft_gays ~ ruralstid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
modelsgroupaffectord [['Transgender individuals']] <- lm(ft_trans ~ ruralstid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)

msummary(modelsgroupaffectord, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

#LGBT rights#
model20 <- lm(lgbtindex ~ ruralstid + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model9)

msummary(model20, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

#rural/small-town ID moderated by place-based identity salience#
#Group based affect#
model21 <- lm(ft_gays ~ ruralstid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model21)

model22 <- lm(ft_trans ~ ruralstid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model22)

modelsgroupaffectordint <- list()
modelsgroupaffectordint [['Gays and lesbians']] <- lm(ft_gays ~ ruralstid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
modelsgroupaffectordint [['Transgender individuals']] <- lm(ft_trans ~ ruralstid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)

msummary(modelsgroupaffectordint, star=FALSE, output='latex', statistic = c("SE = {std.error}", "p = {p.value}"))

#LGBT rights#
model23 <- lm(lgbtindex ~ ruralstid*placeidsalience + place + antistate + contact_LGB + pid7 + ideo7 + White + age + sex + LGBQ + educ + faminc + Evangelical + fattend + region, data = dataframe, weight = weight_post)
summ(model23)

msummary(model23, star=TRUE, output='latex')












